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FOREWORD 


This  is  the  final  report  for  a  three  year  research  effort  entitled 
"Development  and  Implementation  of  Efficient  Sparse  Matrix  Algorithms". 

The  scope  of  this  effort  has  been  very  broad,  covering  the  development  of 
efficient  sparse  matrix  factorization  and  processing  algorithms  and  the 
development  of  local  and  global  optimization  algorithms.  The  topics  form  a 
logical  sequence  in  the  sense  that  the  results  of  each  study  are  supportive 
of  the  successor  studies  in  a  natural  way  for  the  treatment  of  large  sparse 
systems.  Yet,  they  are  autonomous  efforts,  each  having  a  broad  spectrum  of 
applicability  in  its  own  right. 

The  study  on  sparse  matrix  factorization  and  processing  is  reported  in 
"Sparse  Symmetric  Matrix  Processing1  by  P.  S.  Jensen  and  J.  K.  Reid, 

Lockheed  Research  Laboratory  Report  LMSC-D626184  dated  26  May  1978.  This 
report  is  written  as  a  volume  of  three  related  papers  entitled: 

"A  Comparison  of  Two  Sparse  Matrix  Processing  Techniques," 
by  P.  S.  Jensen, 

"A  Package  of  Subroutines  for  Solution  of  Very  Large  Sets 
of  Linear  Finite-Element  Equations",  by  J.  K.  Reid  and 

"A  For*>-an  Virtual  Storage  Simulator  for  Non-Virtual  Computers", 
by  P.  S.  Jensen. 

The  second  paper,  also  available  as  Atomic  Energy  Authority  report  AERE-M  2947, 
documents  the  computer  program  developed  as  a  result  of  this  study.  The 
utilization  of  auxilliary  storage  by  this  program  is  based  on  a  virtual 
memory  concept  that  is  documented  in  the  third  paper.  These  computer 
programs  have  been  distributed  in  the  United  States  and  England  and  requests 
from  Canada  are  being  processed. 

A  preliminary  report  of  the  work  on  local  and  global  optimization  is  pro¬ 
vided  in  "Numerical  Techniques  for  Optimization  and  Nonlinear  Equations"  by 


P.  S.  Jensen,  E.  R.  Hansen  and  H.  M.  C.  Yee,  Lockheed  Research  Report 
LMSC-D630055  dated  29  June  1978.  It  included  a  preliminary  version  of 
"Global  Optimization  Using  Interval  Analysis  -  The  Mul ti -Dimensional  Case" 
by  E.  R.  Hansen,  which  was  subsequently  accepted  for  publication  in  Numerische 
Mathematik.  That  paper  established  the  basis  of  the  computer  program 
documented  in  "GLOBALMIN  -  A  Computer  Program  for  Global  Optimization"  by 
E.  R.  Hansen,  Lockheed  Research  Report  LMSC-D6833Q7 dated  15  November  1979. 
This  program  has  just  been  completed  as  of  this  writing  and  has  not  been 
distributed. 

The  remaining  work  on  local  optimization  is  reported  here.  The  main 
text  covers  the  technical  development  and  test  results.  The  computer  program 
CRATER,  developed  for  the  study,  is  documented  in  the  appendix.  CRATER  is  a 
very  flexible  program  implementing  both  descent  and  quasi-Newton  algorithms 
for  unconstrained,  local  optimization.  The  present  implementation  is 
oriented  toward  small  problems  having  full  Hessians.  An  effort  to  link 
CRATER  to  a  sparse  matrix  processing  system  and  a  non-linear  structural 
analysis  computer  program  (STAGS)  for  large  scale  engineering  analysis  is 
planned  for  the  near  future. 
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ABSTRACT 

A  new  descent  algorithm  based  on  direction  searches  in  the  osculating 
plane  as  opposed  to  the  negative  gradient  is  described.  The  major  practical 
difference  between  the  implementation  of  this  approach  and  quasi-Newton 
algorithms  is  the  use  of  the  Hessian  in  place  of  the  inverse  Hessian.  Thus, 
the  descent  algorithm  is  most  suitable  for  problems  for  which  it  is  practical 
to  recalculate  the  Hessian  a  number  of  times  during  the  solution  process.  A 
number  of  comparative  results  with  popular  quasi-Newton  algorithms  are 
provided.  Extensive  discussion  of  the  implementation  details  is  included 
along  with  documentation  of  a  convenient  computer  program  CRATER  for 
unconstrained  optimization  analysis. 


1. 


INTRODUCTION 


Let  f:Rn  •*  R  be  continuously  differentiable  and  let  g(x)  and  H(x)  denote 
the  gradient  and  Hessian  of  f  at  x.  We  are  interested  in  determining  a 
vector  x*  that  minimizes  f  locally,  i.e.,  such  that  g(x*)  =  0.  For  this 
discussion,  we  consider  only  unconstrained  minimization. 

The  methods  considered  here  are  all  linearization  methods  [13], 
characterized  by  the  use  of  an  affine  approximation 

9 (x)  =  Ak(x  -  xR)  +  g(xk)  (1) 

of  g,  which  is  to  be  valid  for  x  in  a  neighborhood  of  xk-  The 
practical  application  of  this  approximation  is  to  use  its  zero 

=  ?k  -  s!?k> 

as  an  improved  approximation  to  a  zero  of  g. 

We  call  the  change 

~k  =  ~k+l  ‘  ~k 
or 

?k  ‘  -Ak'  9(^k>  (2) 

the  kth  step  vector. 

The  most  famous  linearization  method  is  Newton's  method  which  uses 
Ak  =  ^or  many  Problems,  H  is  very  difficult  or  expensive  to 

-construct.  This  fact  has  given  rise  to  the  widely  studied  quasi -Newton 
methods,  surveyed  in  [3  and  5]  for  example. 

Another  approach  in  determining  Ak  is  motivated  by  considering  x:R  -*■  Rn 
a  continuous,  differentiable  function  of  a  pseudo-time  variable  t  defined  by 


where  x  denotes  the  derivative  with  respect  to  t.  If  x (0)  is  in  a  region 
of  attraction  of  x*,  then  clearly  x(t)  -  x*  as  t  -*  «.  For  a  small  time 
increment  h,  we  may  use  a  Taylor  expansion 

h2  , 

x(t  +  h)  =  x ( t )  +  hx(t)  +  £  x ( t )  +  0(hJ) 

or,  from  (3), 

h2  , 

x(t  +  h)  =  x(t)  -  hg  +  5-  Hg  +  0(lv) 

~  -  -  L  " 

to  approximate  x  by  neglecting  the  high  order  terms  in  h.  Defining 
s  =  x(t  +  h)  -  x(t)  similar  to  (2),  we  obtain 

s  =  ag  +  sHg,  (4) 

? 

where  a  =  -h  and  £  =  h  /2.  We  shall  refer  to  methods  of  this  nature  as 
descent  methods.  In  descent  methods 

A^1  =  al  +  £H(xk) . 

Steepest  descents  is  the  most  well  known  descent  method,  for  which  £  =  0. 

An  important  practical  difference  between  quasi-Newton  and  descent  is 
that  quasi -Newton  uses  H  ^  and  descent  uses  H  to  approximate  s.  If  the 
problem  being  treated  is  such  that  H  is  difficult  to  calculate,  the 
noted  difference  is  of  minor  consequence.  The  appropriate  algorithms  for 
either  model  usually  construct  approximations  to  H~*  or  H  using  low  rank 
(1  or  2)  updates  and  the  two  approximations  are  about  equally  difficult  to 
construct.  If,  however,  it  is  reasonable  to  calculate  H  occasionally 
during  the  solution  process,  (A)  presents  a  clear  advantage  since  it  does 
not  require  a  factorization  of  H.  This  advantage  is  particularly  noticeable 
if  H  is  large  and  sparse.  In  fact,  we  note  in  (4)  that  an  explicit 
'representation  of  H  is  not  required,  i.e.,  only  the  product  Hg  is  needed 
for  various  g.  Thus,  we  can  conveniently  introduce  low  rank  updates  to  the 
current  H  with  minimal  impact  on  an  algorithm  based  on  (4)  by  applying  the 
updates  directly  to  g  rather  than  to  H  prior  to  forming  Hg  (i.e.,  retain 
the  updated  H  in  product  form). 


There  have  been  studies,  e.g.,  [11]  and  [16],  applying  quasi-Newton 
algorithms  to  fairly  large,  sparse  systems.  In  [16],  the  updates  for  the 
inverse  Hessian  were  initially  formed  in  the  conventional  manner  and  then 
zeros  were  inserted  to  conform;  to  the  sparsity  pattern  of  the  problem.  It 
has  apparently  been  shown  to  be  superl inearly  convergent  [5  p  60] .  In  [1], 
the  updates  are  not  directly  applied  to  the  approximate  inverse  Hessian  but 
are  held  and  applied  in  product  form.  After  every  ten  or  so  iterations,  the 
Hessian  is  reformed  and  refactored. 

In  this  paper,  we  consider  the  application  of  a  descent  algorithm  to 
sparse  problems  under  the  assumption  that  it  is  reasonable  to  reform  the 
Hessian  periodically.  As  indicated  in  the  above  discussion,  we  should  not 
expect  as  rapid  convergence  as  with  quasi-Newton  methods  but,  should  expect 
the  average  cost  per  iteration  to  be  substantially  less. 

Ir,  the  next  section,  we  discuss  the  descent  algo'"'thm  used.  It  is  a 
planar  descent  method  in  the  sense  that  a  direction  in  the  [g,  Hg]  plane 
(see  (4))  is  chosen  in  place  of  the  traditional  steepest  descent  direction. 

In  Section  3  we  discuss  the  matrix  updating  schemes  used.  This  topic  has 
been  extensively  discussed  in  the  literature  and  here  we  simply  use  the 
two  approaches  that  have  enjoyed  the  greatest  success  in  the  past. 

Once  the  direction  is  determined,  we  use  a  line  search  to  establish 
the  step  length  as  discussed  in  Section  4.  Thus,  the  function  is  not 
truly  minimized  in  the  plane  [g,  Hg] ,  but  only  along  a  direction  in  that 
plane  determined  on  the  basis  of  (4).  A  study  of  complete  planar  minimization 
at  each  step  remains  to  be  conducted. 

In  Section  5  we  present  some  test  results.  Several  "classical"  problems 
that  have  appeared  in  the  literature  are  presented  along  with  some  new  ones. 

Our  overriding  objective  for  this  effort  was  to  implement  an  effective 
algorithm  for  large,  sparse  problems  for  which  it  is  reasonable  to  construct 
the  Hessian  occasionally  during  a  solution  process.  Although  the  research 
was  directed  toward  descent  methods,  we  considered  it  imperative  that  no 
commitment  to  descent  methods  should  be  reflected  in  the  implementation. 
Furthermore,  because  of  many  variations  in  updating,  step  direction  and  line 
search  that  must  be  considered  in  such  a  study,  we  considered  it  important 
to  use  a  very  flexible  (modular)  design  in  the  implementation  and  include  a 
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convenient  problem  oriented  language  in  order  to  facilitate  a  systematic 
study.  The  resulting  computer  program,  called  CRATER,  is  discussed  in  some 
detail  in  Appendix.  We  feel  that  CRATER  meets  its  design  objectives  very 
well  and  should  prove  very  useful  for  production  development  as  well  as 
additional  studies  in  local  optimization. 


2.  PLANAR  SEARCH 


Instead  of  defining  coefficients  a  and  £  as  in  (4),  we  consider 
predicting  a  and  £  to  minimize  f(x  +  s(a,s)).  Using  a  Taylc'  expansion 
about  x  we  have 

f  ( X  +  s)  :  f(x)  +  sTg(x)  +  ^-STH  (  X  )  s  , 
which  may  be  combined  with  (4)  to  yield 

f(x  +  s)  1  f  +  acQ  +  £a^  +  +  2a£^2  +  —  3^  (5) 

where  we  define 

c.  a  gTH1 g ,  i  =  0,  1,  2,  3. 

The  extrema  occur  at  the  zeros  of  the  partial  derivations  of  (5)  wrt  a  and 
£ ,  which  are  formally  given  by 

,  '  C0C3  (6) 

a  -  2 

°1°3  "  g2 
and  2 

_  °QC2  ~  al 

2 

cla3  "  °2 

2 

if  0^2  f  a 2  .  Otherwise,  we  choose  a1  =  -Og/c^  and  6'  = 

Because  the  approximation  (5)  is  valid  only  for  small  s,  we  have  no 
guarantee  that  the  step 

S  1  =  a  'g  +  6  *Hg  (l7) 

minimizes  f(x  +  s)  over  a  and  £,  or  even  that  f(x  +  s')  <  f(x).  The  following 
theorem  at  least  shows  that  §'  is  in  a  descent  direction. 


Theorem  If  f:Rn  -*•  R  is  twice  continuously  differentiable  with  non-negative 
definite  Hessian  H,  and  at  point  x  e  Rn  the  gradient  g(x)  t  Q,  then  the 
step  s'  defined  by  (7)  is  in  a  descent  direction. 

Proof 

T  2 

We  must  show  that  s'  g  <  0.  If  then  s'  is  parallel  to  the 

negative  gradient,  which  satisfies  our  theorem.  Otherwise  -oz  >  0  by 
the  Cauchy-Schwarz  theorem  and  so  from  (7)  we  evidently  need  to  show 

5  =  (cr3"  c2  )(c5'ao  +  £'°1^  <  0> 

Now  from,  (6)  and  the  definition  of  ; . ,  we  have 

l 

6  =  -(c02a3  -  2cqo1c2  +  a*). 

Def ini ng 

v  -  o gHg  -  c  i  g , 
we  observe  that 

e  =  -v^Hv  <  0.  (8) 

We  have  strict  inequality  in  (8)  because  v  is  not  in  the  null  space  of  H. 
If  it  were,  then  the  denominators  in  the  expressions  for  a'  and  £’  in  (6) 
would  be  zero  (as  would  the  numerators).  That  follows  immediately  from  the 
facts  that  Hv  =  0  implies 

2 

C0C2  C1 

and 

°0  °3  =  °1°2' 

-QED- 

If  we  simply  implement  the  planar  search  algorithm  as  outlined,  and  test 
it  on  the  standard  Rosenbrock  function 
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fR(x)  =  100(x2  -  X]2)2  +  (1  -  x:)2 

starting  at  xT  =  (-1.2,  1),  we  find  the  pleasing  result  that  convergence  is 
achieved  in  only  5  steps.  This  compares  with  3C  steps  for  Powell's  method 
and  18  steps  for  the  Fletcher-Powel 1  method  [6].  Note,  however,  that  this 
comparison  is  not  completely  fair  since  here  we  use  the  true  Hessian.  Also, 
unfortunately,  the  convergence  here  progresses  in  the  rather  erratic  fashion 
shown  in  Figure  1.  It  would  be  somewhat  more  appealing  if  the  function 
decreased  monotoni cally ,  or  nearly  so. 

We  can  arbitrarily  limit  the  factor  by  which  the  function  value  at  any 
point  exceeds  the  value  at  a  previous  point.  For  example,  making  the  factor 
less  than  (or  equal  to)  unity  forces  monotonic  convergence.  In  such  a  case, 
a  multiple  ys'  of  the  step  given  in  (7)  must  be  used,  where  0  <  y  <  1.  Note 
that  a  suitable  y  always  exists  since  s'  is  a  descent  direction  as  snown  in 
the  theorem.  The  value  of  y  must  be  determined  by  a  search  dong  s'.  This 
matter  is  discussed  in  Section  4.  We  shall  see  later  that  imposing  monotonic 
convergence  on  our  planar  search  algorithm  leads  to  a  mildly  slower  conver¬ 
gence  rate  for  the  Rosenbrock  problem. 

If  we  apply  the  basic  method  to  a  quadratic  form 
f(x)  =  xTb  -  \  xTHx 

we  observe  instant  convergence  for  n  =  2  (as  expected)  and  monotonic  con¬ 
vergence  for  larger  n.  For  example,  if  b.  =  n  +  1  -  i  and  H . .  =  i  +  i ,  we 
obtain  convergence  in  5  steps  for  n  =  3  and  convergence  in  15  steps  for  n  =  10. 
The  Newton  method,  of  course  always  gives  instant  convergence  for  this  problem 
anc  quasi -Tewton  methods  theoretically  reouire  no  more  than  n  steps. 
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Figure  1.  Convergence  of  the  Basic  Planar  Search  Algorithm  for  the  Standard 
Rosenbrock  Function  Using  the  Analytic  Hessian. 
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3.  MATRIX  UPDATING 


Techniques  for  updating  the  approximate  Hessian  or  inverse  Hessian  have 
been  very  widely  studied  for  over  25  years.  Dennis  and  More  [5]  provide  a 
recent  survey  of  much  of  the  work  in  this  area  and  Broyden  [3]  provides  a 
very  readable  earlier  survey  of  general  techniques  for  optimization  and  non¬ 
linear  equation  systems.  The  work  by  Broyder,  was  presented  at  the  NSF-CBMS 
Regional  Conference  on  "The  Numerical  Solution  of  Nonlinear  Algebraic  Systems" 
held  at  the  University  of  Pittsburgh  in  1972.  A  highlight  of  that  meeting 
was  a  series  of  lectures  by  Rheinboldt  [13]  that  also  provides  excellent  sur¬ 
vey  information  as  well  as  considerable  theoretical  background. 

For  this  study,  we  used  two  popular  update  procedures,  viz:  DFP  (Davidon- 
Fletcher-Powell )  and  BFGS  (Broyden-Fletcher-Goldfarb-Shanno) .  In  order  to 
present  these  formulas  and  the  form  in  which  they  were  implemented,  we  intro¬ 
duce  some  notation  that  appears  to  be  widely  accepted  among  articles  on  this 
topic. 

As  discussed  in  Section  1,  we  are  seeking  a  point  x*cRn  that  minimizes 
f:Rn  ->  R  locally.  We  use  an  iterative  procedure  which,  given  a  current  point 
xc Rn ,  produces  a  step  SeRn  such  that  x  +  s  is  closer  to  x*  than  x  in  some 
sense.  We  introduce  the  vector 

y  =  g(x  +  s)  -  g(x) 

which  represents  the  change  in  the  gradient  g  of  f  resulting  from  step  s. 
Obviously,  we  would  like  1 |g(x  +  s)'  |  <_  j  g ( x ) '  or  f(x  +  s)  <_  f(x)  or  both. 

Probably  out  of  respect  for  the  extensive  works  of  Broyden  in  this  area, 
the  letter  B  is  usually  used  to  represent  the  approximate  Hessian  and  H  is 
used  for  something  else  (the  inverse  of  B).  In  this  paper,  we  shall  also  use 
B  for  the  approximate  Hessian,  but  H  will  denote  the  true  Hessian  (as  in  Sec.  1). 

If  B  is  the  approximate  Hessian  at  point  x,  we  construct  s,  e.g.,  by 
means  of  the  planar  search  algorithm  discussed  in  Section  2,  and  form  the 
updated  Hessian  using  either  of  the  two  formulas 


Of 


where 


bbfgs  =  6  *  -(JTTT  -  Ty^sT  s  yT 


BDFp  =  B  +  u  zT  +  z  uT 


and 


with 


(y.s)  denotes  the  inner  product  y  s, 

'  =  (y»s)  5* 

v  =  Bs 

u  =  £y  -  v 


>ii  ,  (y,lh 
2{1  (JTiT1 


(9) 


(10) 


It  is  a  simple  algebraic  exercise  to  show  that  (9)  and  (10)  are  formally 
equivalent  to  the  forms  given  on  pp  72,  74  of  Dennis  and  More  [5). 

The  important  features  for  our  application  are:  1_.  The  updates  preserve 
symmetry  (which  is  obvious)  and  2.  BgpgS  s  =  B^pp  s  =  y.  The  latter  feature 
is  basic  to  all  useful  update  formulas. 

For  comparison,  we  also  implemented  quasi-Newton  algorithms  using  the 
BFGS  and  DFP  counterpart  formulas  for  the  inverse  of  B.  Letting  B'1  be  the 
approximate  inverse  Hessian  at  x  we  have  from  [6]. 


where 


bofp  =  B_1  4  TyTiT  -?T  ‘  TJT»T  -  "T 


w  =  B’^y. 


From  [  6]  »  we  have 


where 


BBFGS  =  BDFP  +  -~ 


~  '  (y,§T  ~  ‘  (y.w)  -• 
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(11) 


(12) 


4.  LINE  SEARCH 

Two  fundamental  problems  that  must  be  addressed  in  iterative  optimization 
algorithms  are  selecting  a  step  direction,  as  discussed  in  Section  2,  and  a 
step  length.  Intuitively,  a  length  that  minimizes  the  function  in  the 
direction  of  the  step  seems  most  appropriate.  This  choice  of  length  is 
called  "perfect  iteration"  [3].  However,  that  length  is  fairly  costly  to 
determine  and,  as  we  saw  in  the  example  of  Section  2,  is  not  always  the  best 
choice. 

A  popular  algorithm  due  to  Davidon  [4]  uses  a  cubic  interpolation  of  the 
function  and  its  partial  derivative  in  the  step  direction  correspond!' ng  to 
lengths  of  0  and  1.  If  we  use  length  h  instead  of  1  and  let  f  ,  f^  f ^ , 
fj  denote  the  values  of  the  function  and  its  derivatives  '.along  s)  at  points 
0  and  h  on  the  step  vector  s,  then  cubic  interpolation  suggests  that  the 
length  ;  that  minimizes  f  along  s  is  given  by 

i  =  •  /rii  -  3fg  \  -  1'-i)/3-.2  (13) 

where  ,  ,  ? 

ni  =  (3(fi  _  V  •  h(fi  + 

and  1  1  *3 

r,2  =  >(2(f1  -  f0)  -  h(fj  +  f0  ))/hJ. 

If  f(x  +  as)  is  less  than  f(x)  and  f(x  +  hs),  then  z  is  the  accepted  length 
without  further  searching.  Otherwise  the  process  is  repeated  with  smaller  h, 
e.g.,  |h  or  z.  Since  the  algorithms  usually  used  to  generate  the  step  s 
take  its  length  into  account,  an  initial  h  =  1  is  suitable.  However, 

Fletcher  and  Powell  [6]  suggest  using  conservative  extrapolation 

h  =  min(I,2(f-f*)/f 1 ) 

if  the  value  f*  =  f(x*)  is  known  a-priori . 
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2 

If  f(x)  is  such  that  it  is  practical  computationally  to  use  ;  jgjj  in 
the  place  of  f  in  the  above  algorithm,  some  benefit  accrues  from  the  fact 
that  the  converged  value  is  zero,  viz:  there  is  less  cancellation  in  the 
calculation  of  ' ^  and  and  the  f*  is  always  known  a-priori. 

A  quadratic  line  search  is  also  frequently  used.  Since  four  values 
fQ,  f o ' ,  fj  and  f ^ 1  are  available  for  the  interpolation,  one  can  be  discarded 
or  a  least  squares  fit  can  be  used'.  We  have  chosen  a  quadratic  interpolation 
of  f1  values  that  continues  until  a  zero  (within  a  prescribed  tolerance)  is 
found.  Note  that  this  is  cubic  in  f.  The  processes  are  started  with  f',  f , 

i  ol 

and  f.  ,  where  f comes  from  the  Davidon  cubic  interpolation.  The  amount  of 

Xj 

work  done  in  the  line  search  is  controlled  by  the  tolerance. 


NUMERICAL  RESULTS 


There  have  been  a  variety  of  comparative  results  published  (see  [1.2, 
7,  15]  for  examples)  on  the  unconstrained  optimization  problem.  Several 

problems  appear  repeatedly  in  such  studies  and  seem  to  have  become  defacto 
standard  test  functions.  Occasionally  the  author  of  such  a  study  makes  a 
definite  conclusion  such  as  “The  Fletcher  algorithm  was  clearly  suDerior  to 
all  the  others,  followed  by  the  Davidon-Fletcher-Powell  .  .  . "  in  [7  ]. 

This  author  appreciates  the  value  of  the  test  results  presented  but  cautions 
against  drawing  a  sweeping  conclusion  from  them.  The  problems  involved 
relatively  few  unknowns  (15,  10  and  the  rest  less  than  6)  but  certainly 
represented  a  substantial  complexity.  Questions  relating  to  sparsity  and 
the  use  of  analytic  Hessians  (always  or  occasionally)  were  not  considered. 

Unfortunately,  the  test  results  presented  here  are  also  inadeouate  to 
draw  any  sweeping  conclusions.  They  are  only  intended  to  provide  a 
minimal  indication  of  the  functioning  of  the  algorithms  discussed.  A  more 
comprehensive  study  on  the  behavior  of  these  algorithms  for  the  optimization 
of  structural  panels  will  be  forthcoming.  The  problems  treated  here  are 
described  below. 

5.1  PROBLEMS 

Problem  1.  Rosenbrock  [14] 

A  two  dimensional  problem  that  presents  a  considerable  challenge  for 
numerical  optimization  is  given 

r^(x)  =  olxj2  -  x2)2  +  (1  -  Xj)2 
with  a  =  100  and  starting  at 
xqT  =  (-1.2,  1). 
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In  appearance,  this  function  is  like  a  deep  canyon  with  a  curved,  gently 
sloping  bottom  as  illustrated  in  Figures  1  and  2  (using  a  a  =  10).  Rosenbrocks 
function  converges  to  zero  at  x*T  =  (1,1). 


>  \  /  /./  )  >  ’  n 

,  >  N ;  \  >  >  > 

,  >  %'  ■  &  -  >  ■  • 

>  >  \  ’  ' 

,  >>  >  >  V  '  '  ' 

,■  ,  .  >  \  > 

'  >  >  ,  >  >  (  ,  *  J  '  )  I 

%  % 
''' ' 

*  ’2  ''k  ’>  ' 


Figure  2.  Rosenbrock  Function 


Problem  2.  Wood  [15] 

A  generalization  of  the  Rosenbrock  problem  to  four  dimensions  is  given  by 
f(x)  =  r100(xl’  x2]  +  r90(x3’  x4}  +  19-8  <x2  -  1 ^ x4  '  D 
+  10.1  ((x0  -  l)2  +  (x,  -  l)2) 


7""~r  MB 


» 


i 

i 

I 

? 

I 

i 

i 


where  r  (x)  is  the  Rosenbrock  function  and  the  starting  point  is 
xj  =  -(3,  1,  3,  1). 

Woods  function  converges  to  zero  at  x*T  =  (1,  l,  1,  l). 


Problem  3.  Powell  [24] 

A  problem  that  has  a  singular  Hessian  at  the  point  x*  of  convergence  and 
has  an  extremely  gentle  slope  in  one  direction  near  x*  is  given  by 

f(x)  =  (Xj  +  10x2)2  +  5(x3  -  x4)2  +  (x2  -  2x3)4  +  10(x1  -  x4T 

starting  at 

xj  =  (3,  -1,  0,  1). 

It  converges  to  zero  at  the  origin.  Within  a  distance  e  of  the  origin,  .  ne 

3 

slope  is  0(c  )  along  x^  =  -lOx^,  x3  =  x^. 


Problem  4.  Rayleigh  Quotient 

A  problem  that  does  not  appear  to  have  been  used  much  in  the  past  for 
studies  of  this  sort  is  a  Rayleigh  quotient,  such  as  that  given  by 


f  ( X )  =  ^  riX^Nx 


where 


and 


n  e  n(x)  =  2/x  x 


N  =  diag  (1,  2 . n). 

The  solution  f(x*)  =  1  occurs  at  x*T  =  (1,  0, 
Hessian  of  this  function  are  given  by 


.  . ,  0) .  The  gradient  and 


g(x)  =  n  Nx 
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and 


H ( x )  =  n(N  -  gxT  -  xgT) 

where  N  is  the  shifted  coefficient  matrix  given  by 
N  5  N  -  fl. 

For  application  of  the  planar  search  minimization  algorithm  to  this 
problem,  we  are  concerned  with  the  nature  of  step  vectors  of  the  form  (4). 
Note  that 

T  /  T,.  -Tv 
x  g  =  n(x  Nx  -  fx  x) 

=  2f  -  2f 

=  0 


so  that 

Hg  =  n(Ng  -  gTgx). 

Since  these  results  imply  that 

T  T-  3  T_3 

g  Hg  =  ng  Ng  =  n  x  N  x 

we  see  that  the  vectors  g  and  Hg  are  linearly  independent  unless  x  is  an 
eigenvector,  in  which  case  g  =  0. 

Combining  (4)  with  the  above,  we  obtain  the  expression 

s  =  n(8nN  ^  +  aN  -  sg"^gl)x 

=  (y2n2  +  YjN  +  YqI)x, 

where  Yq  =  n ( 6 ( nf 2  -  gTg)  -  af)  =  nf ( 3£f  -  a)  -  Sn'VN^x 

y ^  =  n(a  -  2snf ) 

-  2« 

Yo  "  n  e- 

For  arbitrary  i  and  j  (i  ^  j)  we  may  write  x  =  x  +  x.e.  +  x.e.,  where  e.  ic 

J  '  J  A  **  ^ 

the  ith  column  of  the  identity  matrix.  Note  that  this  implies  x.  =  x.  =  0. 

«  J 

Note  also  that  the  ith  and  jth  components  of  Nx  are  also  zero.  Letting 
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we  have 


Consequently,  there  exists  a  sequence  of  a's  and  g's  such  that  the  problem 
can  be  solved  in  no  more  than  n/2  steps.  Unfortunately,  the  current  planar 
search  algorithm  does  not  produce  this  sequence  and  somewhat  more  than  n/2 
steps  are  taken. 

The  three  dimensional  Rayleigh  quotient  may  be  graphically  illustrated 
by  parameterizing  it  in  terms  of  coordinates  in  a  plane  passing  through  the 
three  unit  eigenvectors  of  N.  The  result  in  Figure  3,  shows  the  three 
stationery  points  at  (0,0),  (0,1)  and  (1,0). 
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Figure  3.  Illustration  of  a  Rayleigh  Quotient. 


Problem  5.  Quadratic  form 

A  very  simple  quadratic  form  test  problem  was  also  included  in  t1'?  test 
series.  The  form  of  this  was 

f(x)  =  >J(b  +  ^  Hx) 

where 

bi  =  n  +  1  -  i 
and  H  is  diagonal  with 

H..  =  21,  i  *  1,  ...»  n. 

The  solution  f(x*)  =  -  j  £  occurs  at  x*T  =  (n/2,  (n-l)/4 

(n-2  )/6 . l/2n).  1 


5.2  TEST  RESULTS 


Each  of  the  5  problems  described  was  solved  using  each  of  the  3  methods 
(planar  search,  DFP  and  BFGS),  utilizing  each  of  the  3  line  search  techniques 
(cubic,  extrapolated  cubic  and  quadratic).  For  the  planar  search  studies 
each  of  the  3  Hessian  updates  (DFP,  BFGS  and  analytic)  were  used,  making  a 
total  of  75  tests.  The  results  of  these  tests  are  summarized  in  Tables  1 
and  2. 

Not  all  of  the  solutions  converged  within  the  maximum  number  of  ste's 
allowed.  These  cases  are  noted  by  an  *  in  Tables  1  and  2.  The  values  of 
the  functions  and  their  gradient  norms  at  the  end  of  the  iteration  are 
given  in  Tables  3  and  4.  Zero  entries  in  these  tables  correspond  to  con¬ 
verged  results. 
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Table  2 


STATISTICS  for  COMPLETE  QUASI-NEWTON  solutions 
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92 
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i 

ee 

2 
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T 
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2 

73 

74 

2 
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1 
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i 

94 

I 
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A 
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i 

33 

A 

51 

l 

5  i 

I 
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i 
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i 
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1 
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T 
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i 

31 

I 

33 

I 

33 

I 
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I 
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T 

74 

I 
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1  36 

2 
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i 

27 

I 

26 

T 

si 

I 
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Table  3 
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I  I 
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I 

c  UNC 
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I 
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I  ENT 

I 

fun: 

7  ICN 

I 

GRAD] ENT 

I 

fun: 
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I 

CUBIC 
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1 
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I 

3 . 02 

E-06 

1 
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I 

1 

.53 
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l 
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1 
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E-OB 

I 
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I 
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1 
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1 

1 
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I 

6 
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1 
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I 
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T 
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A 
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2 
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7 
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E  00 

I 
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I 

1  .  03  E-03 

i 

0 

i 

* 

4 

A  * 

Table  4 

I — Qli  AS !  -NEwTON  — - - — ------j 

j - DF  p - 1 -  BFGS - 1 

WOOD  -  303  I  FUNCTION  1  GRADIENT  I  FUNCTION  1  GRADIENT  1 

- 1 - 1 - ; - 1 - 1 

CUBIC  I  4.08  E-03  1  2.21  E  00  1  01  01 

EXCUBIC  I  2.75  E-1 1  1  1.15  E-04  2  02  01 

OUAD  I  01  02  01  01 

- 1 - ! - : - 1 - 1 
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6.  CONCLUSIONS 


We  have  presented  a  new  descent  method  that  selects  its  search  directions 
from  the  osculating  planes  instead  of  the  traditional  gradients.  We  have 
described  and  tested  a  variety  of  implementation  considerations  and  made 
extensive  comparisons  with  two  popular  quasi-Newton  algorithms.  The  major 
difference  in  implementation  between  the  descent  and  the  Newton  'or  quasi- 
Newton)  algorithms  is  the  use  of  a  representation  of  a  Hessian  in  the  former 
and  its  inverse  in  the  latter.  When  the  nature  of  a  problem  is  such  that  it 
is  practical  to  work  with  the  inverse  Hessian,  there  appears  to  be  little 
doubt  that  a  good  Newton  or  quasi-Newton  algorithm  will  outperform  a  descent 
algorithm.  When  this  is  not  the  case;  however,  this  new  descent  algorithm 
provides  a  viable  alternative. 

For  both  types  of  algorithm,  we  find  that  the  expenditure  of  a  moderate 
amount  of  effort  in  the  line  search,  such  as  through  the  use  of  our  quadratic 
search  algorithm,  has  a  very  favorable  effect  on  both  convergence  rate  and 
robustness.  We  recommend  somewhat  more  effort  in  this  respect  than  the 
simple  one  or  two  step  cubic  interpolation  that  is  frequently  suggested. 

Two  studies  relating  to  this  work  should  be  made  in  the  future.  The 
first  should  investigate  the  possibility  of  doing  a  plane  search  for  the 
step  length  and  direction  simultaneously  as  opposed  to  obtaining  the  direction 
from  the  plane  and  doing  a  line  search  for  the  length.  The  second  should 
test  the  behavior  of  the  descent  algorithm  on  large,  sparse  problems  and  com¬ 
pare  it  with  that  of  a  good  quasi-Newton  algorithm. 
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APPENDIX 

CRATER  -  AN  INTERACTIVE  PROGRAM 
FOR  OPTIMIZATION  STUDIES 


I 


1.  INTRODUCTION 


CRATER  is  an  interactive  computer  program  designed  for  the  convenient 
solution  of  optimization  problems  using  various  algorithms.  Convenience  of 
operation  is  achieved  by  means  of  a  problem  oriented  language  that  provides 
prompting  whenever  a  user  needs  it  but  does  not  burden  the  experienced  user 
with  unnecessary  questions.  Provision  is  also  made  for  inspecting  inter¬ 
mediate  results  and  changing  various  control  parameters  such  as  the  print 
control  during  an  execution.  The  user  can  even  stop  a  run  and  re-ini tial i ze 
it  at  a  different  starting  point  if  so  desired.  If  a  user  should  become 
confused  at  some  point  in  an  execution,  he  can  type  HELP  for  assistance  or 
STOP  to  quit.  One  or  several  commands  may  be  given  on  one  line.  Finally, 

CRATER  is  forgiving.  A  user  needs  only  to  get  the  first  four  letters  of  each 
command  spelled  correctly.  If  he  fails  at  that,  CRATER  will  politely  ask 
him  to  repeat. 

Convenience  of  problem  setup  is  achieved  through  program  modularity. 

Each  routine  used  by  CRATER  is  designed  to  serve  a  specific  and  rather 
isolated  prupose,  and  is  well  documented  internally.  Three  specific  routines 
with  which  a  user  is  particularly  concerned  are:  USRFNL,  USRGRD  and  USRKES, 
which  define  the  functional,  gradient  and  (optionally)  Hessian  on  which  the 
optimization  procedure  is  to  be  applied.  The  starting  point  can  be  set  at 
zero  or  at  a  variety  of  random  points  using  internal  options  or  it  can  be  keyed 
or  read  in  (free  field) . 

For  transportability,  the  code  is  essentially  written  in  standard  FORTRAN 
66.  We  say  "essentially"  because,  in  fact,  it  is  written  in  a  special  master 
source  code  (MSC)  form  that  includes  directives  for  special,  machine  dependent, 
characteristics  along  with  the  FORTRAN  code  and  resides  in  a  library  (called 
EASY)  developed  at  Lockheed  Missiles  and  Space  Co.  [8].  EASY  is  maintained 
by  a  reasonably  sophisticated  librarian  program  that,  among  other  things,  is 
capable  of  interpreting  the  special  MSC  directives  and  producing  source  code 
that  is  immediately  operational  on  one  of  several  specific  computing  environments. 
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Another  feature  of  MSC  is  that  it  includes  a  structured  program 
documentation  system  that  imposes  a  significant  amount  of  discipline  upon 
developers.  The  librarian  checks  the  documentation  supplied  with  MSC  codes 
and  complains  if  it  feels  that  the  documentation  is  inadequate.  It  also 
extracts  and  tabulates  a  copy  of  the  documentation  in  order  to  facilitate 
library  searches  and  the  construction  of  program  documents  such  as  this. 
Thus,  much  of  the  documentation  appearing  in  the  subsequent  Sections  came 
directly  from  the  CRATER  program  itself. 
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2.  PROBLEM  ORIENTED  LANGUAGE 

The  commands  that  a  user  imposes  upon  CRATER  are  organized  in  a  simple 
tree  structure  as  illustrated  in  Figure  1. 


Figure  1.  Command  Tree  Illustration 


At  level  1,  CRATER  will  list  the  command  options  Cll,  C12,  .  .  .  from  which 
the  user  should  select  one,  C13  say.  CRATER  then  lists  the  level  2  command 
options  under  the  level  1  command  selected,  e.g.,  C24  and  C25,  from  which  the 
user  should  select  one,  and  so  forth.  When  the  user  makes  a  selection  at  the 
end  of  a  branch  in  the  tree,  the  effect  is  recorded  in  the  state  variables 
of  CRATER.  Then  CRATER  moves  back  to  the  next  higher  level  and  repeats  the 
command  options  at  that  level.  The  user  may  then  select  another  command  and 
follow  its  branch  to  the  end.  If  the  user  does  not  wish  to  invoke  any  command 
-options  at  the  current  position  in  the  tree,  he  may  enter  BACKUP  to  move  to 
the  next  higher  level  in  the  command  tree. 

After  some  practice,  a  user  will  tire  of  all  the  prompting  (listing  of 
command  options).  To  avoid  it,  he  may  simply  enter  a  string  of  commands  on 
one  or  more  lines  using  the  line  continuation  flag  $  at  the  end  of  each  line. 
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For  example,  a  user  could  invoke  every  command  illustrated  in  Figure  1  by 
entering  the  following 


Cl  1 ,  C21 ,  C22,  C31 ,  C32,  BACK,  BACK,  $ 
C12 ,  C23 ,  $ 

C13,  C24 ,  C33,  C34,  BACK,  C25,  C35,  STOP 


with  only  a  single  level  1  prompting  of  Cl 1 »  C12,  and  C13.  Note  that  a 
command  may  be  repeated  at  various  times  as  long  as  the  position  in  the  tree, 
at  the  time  the  command  is  issued,  is  appropriate. 

There  are  three  "universal"  commands  that  may  be  invoked  at  any  point 
in  the  tree,  viz:  BACKUP,  STOP  and  HELP.  We  have  discussed  the  BACKUP 
(nr  BACK  for  short)  command  above.  The  STOP  command  simply  terminates  a  run 
and  HELP  will  provide  helpful  instructions  relating  to  the  operation  of  CRATER. 

We  define  two  classes  of  command,  viz:  transitional  and  terminal. 
Transitional  commands  have  "offspring"  commands  at  the  next  level  of  the  tree 
structure  where-as  terminal  commands  do  not.  Thus,  the  terminal  commands 
reside  at  the  ends  of  branches  of  the  tree  structure.  A  terminal  command 
may  require  numerical  data  to  follow  it  but  a  transitional  command  can  only 
be  followed  by  other  commands  in  the  command  stream. 

The  actual  words  used  as  commands  are  simply  established  in  a  data  list 
and  can  be  readily  changed  to  suit  a  users  taste.  See  arrays  COMMND  and 
CMNDPT  in  block  data  deck  L0P91  of  the  code  for  CRATER.  The  itri  word, 
CMNDPT(i),  in  the  pointer  list  points  to  the  i^h  command  word  in  list 
COMMND. 


The  command  tree  structure  is  also  specified  by  data  (see  array  TREE)  in 
L0P91.  However,  the  organization  of  TREE  is  closely  related  to  the  organiza¬ 
tion  of  the  code  itself  and,  consequently,  more  care  must  be  exercised  in 
modifying  TREE.  Users  who  are  quite  familiar  with  the  code  can,  in  a 
straight  forward  fashion,  readily  add  or  delete  branches  on  the  tree.  A 
triple  (0,  S,  W)  is  associated  with  each  node  of  the  tree  where 
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0  is  the  "oldest  offspring"  node  at  the  next  level 


jthe  next  younger  "sibling"  node  at  the  same  level  if  positive  or 
^  1S  (the  parent  node  at  the  previous  level  if  negative 

W  is  the  index  of  the  command  word  associated  with  the  node. 

Up  to  this  point,  we  have  described  the  syntax  of  a  general  POL  organized 
in  a  tree  structure.  The  semantics  of  the  specific  POL  implemented  in  CRATER 
remain  to  be  discussed.  Before  doing  this,  however,  it  is  helpful  to  discuss 
the  organization  of  CRATER  in  order  to  establish  the  setting  in  which  the 
POL  is  couched. 
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3.  PROGRAM  ORGANIZATION 

CRATER  has  a  relatively  small  executive  program  that  oversees  the 
optimization  process.  It  is  supported  by  20  subroutines  that  arry  out 
specific  tasks.  These,  in  turn,  utilize  4  utility  subroutines  frorr  the 
Basic  Linear  Algebra  Library  [10]  and  10  other  utilities  all  of  which  are 
provided  with  CRATER.  The  program,  organization  is  illustrates  in  Table  1. 

Table  1.  Program  Organization. 

Calling  routines  appear  directly  above 
and  to  the  left  of  the  list  of  called 
routines. 


■■■■■ 
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3.1  Purposes  Served  by  CRATER  Routines 

In  this  section  we  provide  a  brief  statement  of  purpose  for  each 
routine  in  CRATER  and  the  utility  routines  supporting  it.  This  section  is 
intended  for  quick  reference  only.  More  documentation  detail  is  provided 
in  Sections  3.2  and  3.3. 


PURPOSE . CRATER 

laboratory  processor  for  STUDYING  local  optimization 

PURPOSE . GTF 

Global  minimization  test  function  evaluation 

PURPOSE . HORNE  R 

Evaluate  a  degree  (n-i)  polynomial  with  coeffs.  c  at  x. 
PURPOSE . LOP05 

accept  AND  INTERPRET  commands  for  local  optimization 
PURPOSE . L0P06 

get  one  user  input  variable  or  command 
PURPOSE . LOP  1  0 

CARRY  OJT  INITIALIZATION  operations  FOR  PROGRAM  CRATER 
PURPOSE . LOP  1  1 

initialize  the  position  vector  for  selected  analysis 

PURPOSE .  ...  LOP  1 2 

INITIALIZE  The  HESSIAN  (UAC03IAN,  matrix 
PURPOSE . L0P20 

CARRY  OUT  ONE  OPTIMIZATION  ITERATION  STEP 
PURPOSE . LOP21 

Evaluate  the  selected  real-valued  functional 
PURPOSE . LCP22 

evaluate  the  selected  gradient  function 
PURPOSE . LOP23 

EVALUATE  the  hessian  OF  the  selected  FUNCTIONAL 
PURPOSE . LOP24 

UPDATE  THE  CURRENT  APPROXIMATE  HESSIAN  MATRIX  (RANK  2) 
PURPOSE . LOp26 

DETERMINE  A  STEP  DIRECTION  S  IN  THE  U‘V  PLANE 
PURPOSE . LOP27 

MOVE  POSITION  X  A  REASONABLE  DISTANCE  ALONG  S 
PURP-SE . LOP30 

CHECK  THE  CONVERGENCE  AND  SAVE  INTERMEDIATE  RESULTS 
PURPOSE . LOP40 

MANAGE  ALL  RESULT  DISPLAY  OPERATIONS  FOR  CRATER 
PURPOSE . L0P81 

MATRIX-VECTOR  MULTIPLY,  U  =  H.V  (H  IS  NORMALLY  A  HESSIAN) 
PURPOSE . LOP91 

PROVIDE  TEXT  DATA  FOR  INTERACTIVE  COMMUNICATION 


PURPOSE . LO  =  92 

PROVIDE  DEFAULT  VALUES  FOR  CONTROi.  VARIABLES  OF  CRATER 
PURPOSE . L0P94 

HOLD  COMMON  DECLARATIONS  GERMANE  TO  MAN-MACHINE  INTERACTION 
PURPOSE . L0  =  95 

HOlD  COMMON  DECLARATIONS  GERMANE  TO  OPERATION  0=  CRATER. 
PURPOSE . L0P98 

DISPLAY  THE  CONTROL  AND  COMMUNICATION  DATA  (PARAMETERS) 
PURPOSE . L0PS9 

GENERAL  PURPOSE  ERROR  HANDLING  ROUTINE 
PURPOSE . USR  c N  L 

CALCULATE  A  SPECIAL  FUNCTIONAL  FOR  USE  by  CRATER 


3.2  Abstracts  of  CRATER  Routines 


In  this  section  we  provide  brief  discussions  on  how  each  routine  in 
CRATER  serves  its  purpose.  Information  on  exactly  how  to  use  each  routine 
is  provided  in  Section  3.3. 
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CPTICS 

”  0  R  ALGOR  IT*-**  7  C  D  E  USED 

9 

6 

PRCA3E 

op  7 : on 

c  OR  PRINTED  jl. T  P  *D%.jVE  (SEE  S  E  w  D  ^  ) 

1  0 

7 

PLCASE 

CP T I  ON 

F'^R  P  L  \j  <  •  z.  j  jU1 

1  1 

e 

LSCASE 

option 

FOR  LINE  SEARCH  TECHNIQUE 

1  2 

9 

X2CASE 

OPTION 

FCR  •  NO’  A  3 S  I  C\ E 0  • 

1  3 

1C 

X3CASE 

OPTION 

FOR  »  NO*  ASSIC.EO  * 

1  4 

1  1 

NV 

PROS wEV 

SIZE  (NO.  I  NOE  RENO  E  N*  '.ARIAS. ES) 

1  5 

12 

MX  I  TR 

VAX.  I* 

era: ION  STEFS  *0  BE  *AKEN 

■ 

:  l 
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1  s 

1  1 

NnOPC7 

WAX.  I" ERA": 

1  7 

14 

L I V  X  1 

LIMIT  PAR4--E 

1  6 

1  5 

L I V  X  2 

LIMIT  PARAMS 

1  9 

16 

nfev 

C..RREN’  NO. 

20 

1  7 

NGEV 

current  no. 

21 

18 

NhEV 

CuRRcNT  no . 

22 

1  9 

NHUP 

C L R R E N *  NO. 

CL 

rc 

N~“  P 

NO. 

2~4 

c  I 

M  *R 

CuRREnT  NC 

25 

22 

HSTAT 

Nj .  OF  06  DEN 

26 

23 

I 2Sr AT 

OP£RAt] ONA . 

27 

24 

I3STAT 

OPERA* : ONAL 

28 

o  r 
_ 

I 4  S  T  AT 

C  R  t  R  A  ^ I DNAl 

rs 

2  3 

I5S7AT 

«  »”  C«M  .  .  V  N  **  w 

30 

27 

iesTAT 

0  -  E  *  A-*  I  SNA. 

3  • 

26 

l5C£Vx 

E .  a \  -  z '  *"Cn 

2  o* 

L  6  C  C  ••  3 

E  „  * *■ •  L  »  \ 

J  J 

DC- 

LB' DVH 

s.av  oo'on 

■5  4 

3 1 

U63GM1 

B^ an-  c:  ‘7 on 

35 

32 

LBCu“»'2 

Blank  common 

36 

33 

LBCCM3 

B.ANK  COMMON 

J7 

34 

L  D  C  0  W  4 

Blank  CCV'C n 

3c 

25 

LE C DM5 

D.  A\w  c:  *  *C\ 

39 

35 

lb  ror-*6 

E-An-  CO  ""ON 

4C 

3" 

LCCCV7 

E.a-.k  CO’  'On 

4  1 

25 

LEOOV5 

B.ANf,  ■  dn 

A  2 

3  3 

L5C0-M9 

E-An*  common 

“  • 

R£; 

.  paramet 

ERS 

43 

1 

ACCENT 

con  v£  p : : nce 

44 

2 

A.=hA 

coefficion-s 

4  "j 

3 

BE  7  A 

T  f-£ 

4£ 

4 

TO  LF AC 

TOLERENCE  :.j 

A  " 

E 

FC3N. 

AN” I C  I  -  *  EL 

A3 

c 

Ri.pr.i3 

•  NO’  ASSIGN 

49 

7 

R  U  PR  7.4 

•  NO*  ASSIGN 

5  G 

8 

RLPRM5 

ACCEPTANCE  C 

5  1 

9 

R  L  P  R  V.6 

pos:*::  .  nor 

C  - 

1  0 

F 

CuRRE'  *  AND 

t  *r 

F  1 

GF  ' ■ 0  F J 

54 

12 

F2 

C  C 

1  2 

GN 

CuRREN'  ano 

56 

1  4 

GN  1 

0-  TnE  GR 

57 

15 

GN2 

58 

16 

GS 

CURRENT  AND 

59 

17 

GS  1 

GRADI ENT- 

60 

IS 

GS2 

61 

19 

SN 

CURREN'  AND 

62 

20 

SN1 

OF  TmE  ST 

63 

21 

SN2 

64 

22 

XI  H 

CURRENT  AND 

65 

23 

XI  HI 

•  \0'  ASS 

66 

24 

X  1  H2 

67 

25 

X2H 

CURR£NT  AND 

69 

26 

X2hi 

•  NOT  ASS 

69 

27 

X2H2 

;  DnS  *3  BE  wSi;  rZk  S~EP  LEN3"r 
: T  £ h  *  NO”  ASSIGNED  * 

[ T E ^  •  NO”  ASSIGNED  » 

C~  FUNCTIONAL  E.A. -A^IOnS 

or  gpa:  : e n t  ;cns 

OF  HESSIAN  e.  -  I  £NS 

C  r  HESSIAN  w ? z A ~  E  S 


STIC  PA F-"E~EF 

s7::  pa-a-.'e^e- 


.  *  n  i  -  : 
‘  r  c  »  c  : 


NC*  A5S 


v.  2  - .  c  :  ^  ;  £  £  3  ~  n  p 


ASSIGNED  • 

-IRE  *  BITS  A- £  FP.EE  , 


- ; t e t : a  fop  : 

(IP  CA.Cw'.L 


»  A  -  w  £  S 


.  ^  -v- 

***'’  "-1  . '  *"j  I 


PREVIOUS  .  ALuES 
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The  INDEX  of  THE  USER  COMMAND  IN  THE  COMMA'. 


■REE  (AT  THE  CURRENT 
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LEVEl)  IS  DETERMINED  :CR  THE  CALLING  PROG"-''.  T  -i  1 5  INDEX  IS  ASS_  ED 
TD  INDICATE  TO  THE  CALLING  PROGRAM  WHAT  IT  IS  EXPECTED  TO  DC. 

ABSTRACT . . . . LOPOS 

A  jSER  MAY  I N p l1  t  SEVERA.  ITEMS  ON  ONE 
IS  STACKED  IN  A  COMMON  ARRAY  /^OPJNP/  AND 
CCF  BY  LOR'06.  ONE  FOR  ETCH  CA.L.  IF  TnE  S 
DA  i — S  FOR  MORE  INPuT, 

ABs’RACT. . . . LCP10 

this  routine  carries  out  initialization  commands  issued  by  a 

USER  IN  ACCORDANCE  w!*H  THE  COMMAND  LANGUAGE  TREE  DF  TnE  CRATER 
OP* I V I D  AT  I  ON  PROGRAM. 


stdase  init:a.iza*:on 

STANDARD  FOR  FDStN BROCK  BANANA 


6  RAND 

ZERO 


AEB'RACT.  .  .  .  LO-M  2 

INITIALIZE  T-i £  HESSIAN  MATRIX  ACCORDING  TO  thE  CR I  *  E  R I  O'.  GIVEN  I 
PARAMETER  IriCASE. 

INCASE  INTERPRETATION 

1  ANA.YTICAl  HESSIAN  !ROl*:nE  LOr2?; 

«  iDtNTITv  MATRIX 

w  R c. AD  F«OM  AuX.  STORAGE 


SENSE;  THAT  evaluate  the  CORRESPONDING  GRADIENTS  AND  HESSIANS.  ThE> 
ARE  CALLED  lOP22  AND  LQR23. 

CURRENTLY.  THE  cJNCtIONAlS  DEFINED  ARE  AS  FD..DWS 

case  funct i on a  l 

1  RDSENSROCK  BANANA  (S'ANDARD) 

2  SIMPLE  ELlIPSDIO 

3  special  (user  supplied) 

a  pd* e l l  Function 

ABSTRACT - LOP22  tlf.t  LK  tUVi  i 

SEE  ABSTRACT  for  LGR21  OJf  1  #*£;.  BO  DW 

AESTRACT. . . . L0P23 

SEE  ABSTRACT  FOR  L0P21 
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ABSTRACT. . . . lOP24 

TriE  FuLl  HESSIAN  OP  INVERSE  HESSIAN  IS  U  =  : 
R'N*.  2  DFP  ( DA V I  DON- F  lETChER- ROWS  L  l )  FOR7U.A  S 

*j . E  . o E \ n I  s ,  u R  •  and  j.  j.r.'CRE,  quasi-ne/.tdn  vet: 
1  9, 1  ( ‘ 977 )  46,83. 


ED  USING  ThE 
SEE.  E .G.  . 

S.  S I A  V  REVIEW. 


IN  PRACTICAL  APPlICA* IG\S  CaITh  5= A=SE  Hi.  n  S~D_.D  NO*  Si 
r.  £  D "  t v  MDDI  FI  ED  .  INS'EAD.  THE  \'ECTCSS  i  AND  -  ’AND  iDA.E^S 


e  . 

S-OJlD  EE  re 

'AINED 

FDR  USE 

IN  S-5SBC - B NT  CALCw.A-IONS  C-* 

H*X  F;p 

arbitrary  X 

AESTRAC 

T. . . . L0P26 

Z  E  ~ 

ERMINE  S  = 

A  L  *  U 

-  3"  »  V 

CDDaRS 

AT  A  POINT  N 

E  A  -  £  ^ 

TD  X-rS 

T  -1  -  \  x.  J  ^£3E\T  i  .  '  E  r  :  ap-^IEE 

T  3  ‘J  s 

G  AND  V  = 

H»  G . 

a  £  A  E  G 

IS  *  h  E  c.  —Cn*  gr-disnt  . 

ABSTRACT. . . . lD  =  2? 

FOR  THE  CUBIC  search,  the  step  length  is  accepted  when  the 
FUNCTION  value  does  NOT  exceed  TOLFAC*FO,  WHERE  TOLFAC  is  a 
PARAMETER  and  FO  IS  the  initial  FUNCTION  VALUE.  for  QUADRATIC 
search,  A  ZERO  OF  The  PARTIAL  DERIVATIVE  of  The  FUNCTION  IN  The 
STEP  DIRECTION  IS  sought  UNDER  THE  CONTROL  OF  PARAMETER  RLPRM5. 


ABSTRACT _ lC?30 

DESIRED  STATISTICS  ARE  RECORDED  AND  A  DON »  C - GE  VC  E  C-CC 
MADE.  THE  CONVERGENCE  CRITERION  IS  ThAT  t-.s  »  hRV  C:  *-E  GRA 
NSC.ICCA3LE  IN  CCVPARISDS  a  I TH  T.-.E  ACCEPTANCE  PARAMETER  AC. 


L0P4C 

CR  AT  ER  CCV'AND  INTERPRETER  *  C  DECI-hER  JiER  CCA-'ANDS 

out  those  commands* 

ASSTRACT. . . . lC=8i 

special  matrix-vec*cr  m.lt:c>.y  rc.’In 
*u:s  routine  takes  into  accoun'  the  s:eci 

C f* I " I  2  AT  ION .  I.E..  H  IS  A  SASIC  IS  =  AR5E> 

C r  A.SCCCIA'lD  RANK-1  l’PDA'E  V  SC'jRS. 

currently,  h  is  uust  a  fu.v  matrix  with  the  .•  pdate  s  exp.icit„v 


ABSTRACT  .... L0P91 

this  ESTABLISHES  ThE  TEXT  FOR  INTERACTIVE  COMMUNICATION  BETWEEN 
THE  OPTIMIZATION  PROGRAM  CRATER  AND  A  USER. 

ABSTRACT. . . . lO°S2 

Trris  ESTABLISHES  All  0^  THE  DEFAULT  VA.uES  F-=.  CONTROL  VARIABLES 
USED  IN  PROGRAM  CRATER. 

ADS'RACT - uCP9e 

Systematically  disp.ays  the  con'ents  cf  the  hE'  lase.lED  common 
BLOCKS. 


ABSTRACT . . . . lC=99 

AN  IN’ERNAu  table  IS  MAINTAINED  WHICH  INDICAtES  whA*  AC'ION  is 
TO  BE  TAKEN  FOR  EACH  OF  THE  L0°  ROUTINES.  THIS  ROUTINE  CARRIES  OUT 
THE  INDICATED  ACTION 


& 

s 


UI3  tA&M  Li  QVALi  l  i  *  . 

rwM  (sun  woLitsttau  so  dlc 
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3.3  Directions  for  Using  CRATER  Routines 

In  this  section  we  specify  the  parameters  and  other  considerations 
involved  in  using  each  routine  in  CRATER  and  the  utility  routines  supporting 
it. 


_  5A  St  . CRATER 

S:V  =  LY  EXECUTE  CRA'ER  AND  INTERACTIVELY  S.PF.Y  COMMA :.CS  AND 
DA'A  AS  REQUESTED.  THE  TREE  STRUCTURE  OF  *HE  CCMMAND  lANO.ACE 
IS  EX=.A\EC  IN  RESPONSE  TO  THE  HElP  COMMAND. 

Tr.E  COMMANDS  APPROPRIATE  FCR  'HE  CURRENT  PQSI'ICn  IN  *-iE  CCV'-'.C 
TREE  ARE  DISPLAYED  AND  ThE  USER  HAS  THE  FC. LOwING  CD'ICNS 

1.  T'.PE  ONE  OF  THE  COMMANDS. 

2.  * > PE  THE  WORD  HELP  FOR  GUIDANCE- 

3.  -YPE  The  WORD  EAC-v-A  TO  3i’  TO  PPC\I0jS  POSITION  IN  ThE  'REE. 

4 .  TYPE  ThE  WORD  STC=  TO  TERM I NAT £  ThE  ExSC.'ION. 

IN  ORDER  TC  DC  LOCA.  I  N I V 1 2  AT  I  ON  ON  A  USER  SUD  =  liED  FjnC'IOna.. 
'HE  .SER  Ml  ST  Supply  SUBROUTINES  ,.  T-E  cO. LOWING  CAlLING  SEQUENCES 

CALl  USR-N-l N , c . X .  W) 

CALL  U  S  R  G  R  0  {  N  ,  3  .  X  ,  *< ) 

CALL  USRHES(N.H.X.W) 

WHERE  F  IS  THE  SCALAR  FUNCTIONAL  VALUE  AT. POSITION  X 

G  IS  THE  VECTOR  GRADIENT  AT  POSITION  X  : LENGTH  Ni 
H  IS  THE  N  BY  N  HESSIAN  HA'R’X  a*  POSITION  x 

n  is  the  no.  independent  variables  in  the  Functional 

X  IS  THE  CURRENT  POSITION’  VECTOR  (INDEPENDENT  VARIABLE) 

W  IS  THE  WORKSPACE  (SCRATCH)  VECTOR  (LENGTH  N) 

The  DETAILS  OF  THE  PROGRAM  FLOW  ARE  ALMOST  INTIRELY  CONTROL. ED  BY 
-ES  ASSIGNED  TO  PARAMETERS  IN  LA5E -  LED  COMMON  'SEE  ASS'RACT). 

ft  may  be  changed  3y  specify ing  'he  type,  index  ant  .a.^e,  e.c-  . 
E3EP.2.3  tc  SET  THE  SECOND  INTEGER  TC  T-E  VA.lC  3.  This  TYPE  OF 
K.  IS  DONE  AT  THE  PARAMETER  NODE  OF  THE  COMMAND  TREE- 

CO’MON  VARIABLES  AND  arrays  ARE  ESTABLISHED  IN  EXTERNAL  PROCED¬ 
URES!  COMDECKS  OP  INCLUDE  DECKS)  lOP94  AND  -0?95,  AND  IN  EASV  PROCED¬ 
URE  MTNQT.  THESE  PROCEDURES  MUST  BE  MADE  AVAILABLE  FCR  T nE  COMP¬ 
ILATION  OF  THE  REST  OF  THE  CRATER  ROUTINES. 

THE  amount  OF  PRINTOUT  IS  CONTRO.LEO  by  A  PR1Nt  LEVEL  ES'AElI SHED 
BY  MEANS  OF  A  PRINT  COMMAND.  EACH  LEVEL  (C  TO  101  INCLUDES  THE  CUTPUT 
OF  LOWER  LEVELS.  A  ROlSH  DESCRIPTION  OF  ThE  OU'p-T  INTRODUCED  AT 
EACH  LEVEl  IS  GIVEN  BE’.Cw.  THE  NUMBERS  IN  PARENTHESES  INDICATE 
FROM  WHICH  ROUTINE  THE  OUTPUT  COMES. 

C.  SUMMARY  data  (40)  AND  FATAL  ERRORS  (59’ 

1.  final  POSITION  (SOLUTION)  (40),  NOTE  NONDESCENT  STEFS  (27) 

2.  CONVERGENCE  HISTORY  (40),  NGNCONVERGED  STEP  MULTIPLIERS  (27) 

3.  NOTE  OCCURRENCE  OP  COMPLEX  ROOTS  IN  QUADRATIC  LINE  SEARCH  (27) 

4.  DISPLAY  INITIAl  POSITION  VECTOR  (11) 

5. 

6.  FUNCTIONAL  VALUE  as  CALCULATED  (21).  DEGENERATE  PLANES  (26) 

7.  GRADIENT  (22)  AND  POSITION  (20)  VECTORS  AS  CALCULATED 

8.  STEP  VECTOR  AS  CALCULATED  (26),  ERROR  ACTION  CONTROLS  (98) 

9.  HESSIAN  AS  CALCULATED  (23).  COMMAND  DA'A  (93) 

10.  LINE  SEARCH  DETAIL  (27),  CONVERSATIONAL  DATA  (93) 

USAGE  . LOP05 


-  FUNCTIONAL  evaluation 

-  C-RADIEN*  evaluation 

-  HESSIAN  E.A.UA'ICN  (QPTIONA.) 
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USAGE  . LOP22 

CALL  LQP22 (  N  .  G  .  X  ,  w  ) 


ARG. 

TYPE 

PURPOSE 

N 

INT. 

LENGTH  0=G,  X  AND 

,v  (  ARCE  _  Zf» 

G 

REAl 

R  tSJ  I  N'G  G  R  A  G  *  E  \  • 

VECTOR 

X 

REAl 

INDEPENDENT  VASIAS 

-  E  E  Z  ~  GR 

w 

real 

WORK  SPACE  ( LENGTH 

N) 

<J  S  A  G  £ 

LOb23 

CA 

L  L  L0p23 ( N , H , X . ft ! 

A FG  . 

TYPE 

PUP  =05  E 

N 

H 

IN’  . 

NG. 

\a  r :  a: _e: 

REAL 

N  E  •  \  hicb.A  v 

**  »  A  .  V-  -  *'h  / 

X 

real 

I  ND1  r  E\GZ‘.*  VAR1AE 

_ E  V£ GTC“ 

w 

real 

ft  OR  r,  SPACE  Cc  LEND 

T  ri  N 

USAGE  . LCf'24 

CALL  LOp2 A (  N  ,  H  ,  S  .  G  ,  V.  I 
A RG •  TYPE  PUP  ^CSc 

n  int.  PscE.Ev  size  (no.  indep enoen'  .asi-e.ee 

h  pea.  n  e •  n,  - _ -:.::a’- 

S  REA.  STEP  LSC’CP  i lENSTh  Ni 

G  real  Current  gradient  vec'Gr  iivediate.  rep^aClc  ev  o  - . ,  > 

W  REA.  PREVIOUS  GRADIENT  VECTOR  A\U  ».  oRK  SFA.E  (Sl  0  P*3 1 

USAGE  . LOP26 

CALL  L0~26f  M , S,U,\  , H ) 

ARG.  TYPE  PURPOSE 

N  INT.  PRGc.EV  SICE 

S  REA.  S*  E 1  VECTCR  RROC-CEO  -.ERE  FCR  NEX'  APPLE/.  .ENG*-  N' 

U  real  COMPONENT  of  P.ANE  ON  ft-IC-  '0  SE- P C-  FOR  S  ( LENOT-  !, 

V  REAL  COMPONENT  of  P.ANE  ON  ft-IC-i  ‘0  SEARC-  *Z-  S  i  .ENG*- 

H  REAl  N  3t  N  nEESIAN  V  A  T  R  I  X  (  0  U  R  R  E  N  .  v  r.L  - 

USAGE  . lO=27 

CALL  L0P27 ( N, X , S , G, W ) 

ARG.  TYPE  PURPOSE 

N  INT.  PROBlEV  DIMENSION 

X  REAL  CURRENT  POSITION  (INDEPENDENT  VARIABLE)  (LENGTH  N) 

S  REAL  STEP  DIRECTION  FOR  NEW  POSITION  (LENGTH  N) 

G  REAL  GRADIENT  VECTOR 

W  REAL  WORK  SPACE  (LENGTH  N«2) 


THE  technique  for  The  line  SEARCH  is  CONTROLLED  BY  THE  integer 
parameter  lscase  as  follows: 

1  T  -  CUBIC  (My  FORMULA) 

1  2  -  CUBIC  (my  FORMULA)  WITH  1ST  POINT  EXTRAPOLATION 
LSCASE  «  <  3  -  QUADRATIC 

I  4  -  CUBIC  (DAVIDON  FORMULA)  WITH  1ST  POINT  EXTRAPOLATION 
I  5  -  CUBIC  (DAVIDON  FORMULA) 

THIS  ROUTINE  USED  PARAMETERS  ALPHA,  TOLFAC  AND  RLPRM5  ALSO 
(SEE  THE  ABSTRACT) . 

USAGE  . LCP30 

CALL  LCP301CONV) 

ARG.  TYPE  PURPOSE 

CON  V  LOG.  SET  TRwE  IF  AND  ONLY  IF  THE  PROCESS  hAS  CONVERGED 
AS  REQUIRED 


USAGE  . LOP40 

CALL  LOP40  (NO  ARGUMENTS) 


RESPONDS  TO  DISPLAY  COMMANDS  ISSUED  BY  T ni  -SER 


USAGE  . LCP61 

CALL  LCP31 (N.U.H.V) 

ARG.  TYPE  PURPOSE 

N  INT.  DIMENSION  OF  SYS*  EM 
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4.  COMMAND  INTERPRETATION 


In  this  section  we  discuss  the  specific  commands  implemented  in  the  POL 
(problem  oriented  language)  of  CRATER.  This  discussion  draws  heavily  on  the 
material  presented  in  Section  2  and  is  probably  unintelligible  for  the 
reader  who  has  not  first  read  that  Section.  A  general  familiarity  with 
Section  3  will  be  helpful  for  the  reading  of  this  section  but  is  not  quite 
as  crucial . 

Both  the  tree  structure  and  the  . ommand  words  used  in  CRATER  are  pro¬ 
vided  in  Table  2.  As  mentioned  in  the  introduction,  only  the  first  four 
letters  of  each  command  are  required  from  the  user.  For  each  command  word 
C,  the  "offspring"  command  words  (command  word  options  at  the  next  level 
of  the  tree  that  are  associated  with  C)  appear  indented  and  directly 
under  it.  Table  1  is  constructed  by  CRATER  in  response  to  the  commands 
DISPLAY,  COMMANDS  (or  DISP,  COMM)  issued  at  tree  level  1.  Note  that 
several  command  words  such  as  PRINT,  READ,  DISPLAY,  etc.  are  repeated. 

They  appear  only  once  in  the  command  word  list  COMMND  but  are  associated 
with  several  nodes  of  the  tree. 

CRATER  indicates  the  current  level  in  the  tree  whenever  it  lists  the 
command  options  for  the  user.  For  example, 

?  PLEA  TE  TYPE  ONE  OF  THE  FOLLOWIN'?  COMMANDS  ... 

FUNCTION  PPINT  ADDITION  HE  DTI  AN 

»APPf>iFTFP  EVALUATE  METHOD  DEA»CH 

is  a  list  of  level  2  options  associated  with  INITIALIZE.  If  a  user  has 
lost  track  of  what  options  he  has  selected  up  to  a  certain  point,  he  can 
back  up  to  level  1  (using  one  less  BACKUP  command  than  the  current  level 
-number)  to  get  to  level  1,  and  issue  commands  DISPLAY  OPTIONS.  A  typical 
response  is  as  follows: 


TABLE  2.  Command  Table  for  CRATER 


1 

initialize 

2 

PRINT 

3  . 

function 

4 

.  ROSENBROC 

5  . 

ellipsoid 

6 

SPECIAL 

7  . 

POWELL 

8  . 

WOOD 

9  . 

POSITION 

10 

.  ROSENBROC 

1 1 

.  RANDOM 

12 

ZERO 

13 

.  READ 

14  . 

.  POWELL 

15  . 

WOOD 

16 

.  RAYLEIGH 

17  . 

HESSIAN 

18  . 

INITIAL 

19 

.  .  analyti 

20  . 

IDENTIT 

21 

READ 

22 

UPDATE 

23  . 

.  .  analyti 

24  . 

DFP 

25  . 

.  .  BFuS 

26  . 

method 

27  . 

PLS1 

28  . 

DFP 

29  . 

BFGS 

30  . 

SEARCH 

31  . 

.  CUBIC 

32  . 

.  EXCUBIC 

33  . 

quadratic 

34  . 

parameter 

35  . 

.  INTEGER 

36  . 

REAL 

37  . 

.  DISPLAY 

38  . 

evaluate 

39  . 

ITERATE 

40  . 

DISPLAY 

41 

results 

42 

.  PRINT 

43  . 

.  plot 

44 

.  .  termina 

45 

.  .  SYSTEM 

46 

OPTIONS 

47 

COMMANDS 

48 

STATE 
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current  options 


1 

POSITION 

WOOD 

2 

FUNCTION 

WOOD 

3 

INITIAL 

HESSIAN 

IDENTITY 

4 

UPDATE 

HESSIAN 

DFP 

5 

PRINT 

9 

6 

METHOD 

DFP 

7 

SEARCH 

EXCUSIC 

As  discussed  in  Section  2,  two  classes  of  command  are  used,  viz: 
transitional  and  terminal.  Tne  terminal  commands  reside  at  the  ends  of  the 
branches  and  certain  of  them  require  input  data  immediately  following  them. 
Note  that  any  command  in  Table  2  that  does  not  have  another  command  directly 
beneath  and  to  the  right  (indented)  of  it  is  a  terminal  command.  Thus,  for 
example,  PRINT,  ANALYTIC,  and  ITERATE  are  terminal  commands.  Prompting  is 
provided  for  terminal  commands  requiring  input  data,  such  as  READ. 

4.2  Coma-, and  Semantics 

The  interpretation  of  most  of  the  commands  is  self  evident.  Here  we 
shall  run  through  the  commands  in  Table  1,  pointing  out  special  considerations 
as  we  go. 

4.1.1  Initialization 

The  first  command  to  be  entered  is  usually  INITIALIZE,  however,  a  user 
might  enter 

DISPLAY  OPTIONS  BACK  INITIALIZE 

in  order  to  find  out  what  the  default  options  of  CRATER  are  before 
initial izing.  Except  for  EVALUATE,  which  should  be  the  last  command  used 
during  initialization,  the  order  in  which  the  initialization  commands  are 
imposed  is  flexible.  EVALUATE  establishes  the  function,  gradient  and  Hessian 
.values  at  the  initial  position  and  thus  requires  that  they  be  previously 
defined. 

There  are  three  standard  test  functions  (ROSENBROCK,  ELLIPSOID  and 
POWELL)  built  in  and  a  l’nk  SPECIAL  to  a  user  defined  function.  The  terminal 
commands  ELLIPSOID  and  SPECIAL  require  the  number  of  unknowns  in  the  problem 


to  be  entered  as  data  (free-field  integer).  Since  it  only  makes  sense  to 
select  one  function,  CRATER  generates  a  BACKUP  command  internally  after  a 
selection  is  made  for  the  convenience  of  the  user.  This  innovation  applies 
also  to  POSITION,  INITIAL,  UPDATE,  METHOD,  SEARCH,  and  PLOT  commands. 

The  PRINT  command  sets  the  amount  of  operational  detail  to  be  printed 
out  during  the  analysis.  It  must  be  followed  by  a  print  level  (integer) 
ranging  from  0  to  10.  The  effects  of  the  various  print  levels  are  described 
in  the  CRATER  USAGE  description  in  Section  3.3.  Typically,  the  print  level 
is  set  high  (8  -  10)  for  the  first  couple  iterations  and  then  lowered  to 
(0  -  5)  for  the  rest. 

The  POSITION  command  should  fol low  the  FUNCTION  command.  The  three 
standard  initial  positions  are  built  in  corresponding  to  the  standard,  built 
in  functions.  In  addition,  a  random  or  zero  initial  position  can  be  generated 
internally  and  the  user  has  the  option  of  entering  a  special  initial  position 
after  the  READ  command.  If  the  size  of  the  function  has  not  been  established 
before  entering  POSITION,  it  does  not  know  how  large  the  initial  vector 
snould  be.  A  user  can  obtain  a  variety  of  random  initial  positions  by 
repeating  the  command  sequence  POSITION  RANDOM. 

The  method  for  constructing  the  HESSIAN  must  be  established  both 
initially  and  for  update  purposes  during  the  analysis.  The  Davidon-Fletcher- 
Powell  (DFP)  and  Broyden-Fletcher-Goldfarb-Shanno  (BFGS)  update  formulas  are 
built  in,  both  for  updating  the  Hessian  and  the  inverse  Hessian.  For  this 
study,  an  analytical  initial  Hessian  is  used  with  either  the  DFP  or  BFGS 
update  formula.  Periodically,  during  the  analysis,  the  analytical  Hessian 
is  reconstructed  by  means  of  the  command  sequence 
INITIALIZE,  EVALUATE,  BACK 

and  then  the  iteration  is  continued.  It  is  possible  to  have  that  commanc 
sequence  generated  internally  by  the  iteration  portion  of  the  program  (on 
the  basis  of  a  suitable  decision  heuristic)  but  this  has  not  been  done. 

Most  of  the  operational  "state"  variables  are  maintained  in  a  parameter 
list  described  in  the  CRATER  documentation  of  Section  3.2.  The  PARAMETER 
command  gives  the  user  access  to  these  parameters  both  to  see  what  values  they 
have  (DISPLAY)  and  to  change  them  (INTEGER  and  REAL).  For  example,  to  set 
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the  4th  floating  point  parameter  to  5000,  one  may  enter 

REAL  4.5.E3 

or  seme  FORTRAN  equivalent  representation  of  5000.  It  is  essential  to  under¬ 
stand  the  purposes  of  the  parameters  (see  Sec.  3.2  under  CRATER)  before  doing 
any  modifications. 

The  methods  currently  implemented  are  the  planar  search  with  one 
dimensional  line  search  (PLSI)  descent  method  and  versions  of  the  famous 
Davidon-Fl etcher-Powell  (DFP)  and  Boyden-Fletcher-Gol dfarb-Shanno  (BFGS) 
quasi-Newton  methods.  In  the  current  implementation,  the  Hessian  updates  are 
applied  directly  to  the  approximate  Hessians  rather  than  being  held  in  product 
form.  Consequently,  it  is  not  appropriate  to  use  an  analytical  initial 
Hessian  with  either  the  DFP  or  the  BFGS  algorithms  since  the  requisite  matrix 
inversion  algorithm  is  not  provided. 

A  choice  of  the  Davidon  cubic  interpolation  line  search  (CUBIC)  without 
extrapolation  or  with  it  (EXCUBIC),  or  a  quadratic  line  search  (QUADRATIC) 
is  provided.  If  EXCUBIC  is  selected,  then  the  expected  value  of  the  converged 
function  must  be  supplied  as  data. 

4.1.2  Iteration 

The  iteration  process  is  trivial  to  invoke  in  comparison  with  initiali¬ 
zation.  The  user  simply  enters  ITERATE  followed  by  the  maximum  number  of 
iteration  steps  to  be  taken.  As  mentioned  earlier,  it  is  convenient  to 
initially  set  the  print  level  high  and  iterate  only  2  or  3  times  in  order  to 
see  that  the  process  is  going  properly.  The  user  can  then  return  to 
initialize  and  lower  the  print  level  if  everything  is  ok,  or  else  change 
the  initial  conditions  in  order  to  improve  things. 

4.1.3  Display 

Besides  intermediate  results  that  can  be  obtained  during  the  iteration 
process  by  setting  the  appropriate  print  level,  CRATER  provides  a  variety  of 
summary  results  via  the  DISPLAY  command.  Values  for  the  function,  the 
gradient  norm,  the  step  norm  and  inner  product  of  each  new  gradient  with  the 
previous  step  (0  for  perfect  iteration)  are  tabulated  during  iteration  and  can 
be  printed  or  plotted  (on  Univac  computers)  in  the  display  mode.  The  plotting 
system  is  machine  dependent  and  would  undoubtedly  have  to  be  modified  for 
different  computing  environments.  Terminal  plotting,  of  course,  is  only 
applicable  if  the  user  has  a  graphics  terminal. 
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In  addition  to  computational  results,  select  options  chosen  during 
initialization  and  the  complete  command  table  (see  Table  1)  can  also  be 
generated.  Should  something  go  wrong  with  the  computation,  the  program 
state  variables  may  also  be  displayed  for  debugging  purposes.  This  output 
tends  to  be  rather  long  and  is  seldom  needed. 


4.2  Batch  Processing 

The  problem  oriented  language  of  CRATER  provides  a  convenient  means 
for  setting  up  batch  jobs  as  well  as  interactive  ones.  The  following  is 
typical  of  such  a  job  in  a  fairly  readable  form. 


HELP  $ 

INITIALIZE  $ 

PRINT  10  $ 

FUNCTION  ROSENBROCK  $ 
POSITION  ROSENBROCK  $ 
HESSIAN  $ 

INITIAL  ANALYTIC  S 
UPOATE  ANALYTIC  $ 
BACK  $ 

METHOD  PlS!  $ 

PARAMETER  $ 

REAL  4.5000.  $ 

BACK  S 

search  cubic  $ 

EVALUATE  $ 

parameter  s 

DISPLAY  $ 

BACK  $ 

BACK  $ 

DISPLAY  $ 

OPTIONS  $ 

BACK  S 

ITERATE  1  S 

INITIALIZE  S 

PRINT  5  $ 

BACK  $ 

ITERATE  200  $ 

DISPLAY  $ 

RESULTS  S 

PRINT  S 

STOP 


WE  FIRST  ILLUSTRATE  the  HELP  INST  R'JC*  I  ON  • 

this  is  a  typical  problem  definition  in  the 

CRATER  Pol  (PROBLEM  ORIENTED  lANGJAGE).  THE 
PRINT  LEVEL  IS  SET  HIGH  TO  OBTAIN  ALL  OUT¬ 
PUT  DURING  INITIALIZATION.  THE  STANDARD 
ROSENBROCK  FUN.  IS  TO  BE  SOLVED  USING  rs 
STANDARD  initial  POSITION  VECTOR.  THE  IMT- 
IA^  AND’  UPDATE  HE  SS I  AN'S  ARE  TO  BE  ANAlytIC, 
demonstrating  the  BASIC  algorithm  discussed 

IN  TEXT.  THE  PlANAR  SEARCH  1  (PlSI)  ALGORITHM 
IS  TC  BE  USED.  REAL  PARAMETER  4  (TOLFAC>  IS 
SET  ’0  5000.  IN  ORDER  to  free  all  MOND'CN- 
i c i t  »  constraints,  we  shall  use  the  cubic 
sea=:h  in  theis  setting. 

after  the  Initial  function,  gradient 

and  HESSIAN  EVALUATION.  WE  DISPLAY  the  PAR¬ 
AMETERS  FOR  DOjBlE  CHECKING  LATER.  NOTE 
That  we  always  enter  Back  to  return  TO  the 
PREVIOUS  tree  level  (INDENTATION  LEVEL).  WE 
choose  to  obtain  a  summary  of  the  options 

SElECTEO  On  the  OUTPUT  for  HANDY  REFERENCE. 
NEXT  WE  CARRY  OUT  ONE  ITERATION  WITH  LOTS 
OF  PRINTOUT  AND  THEN  GO  BACK  TO  INITIALIZE 
IN  ORDER  TO  REDUCE  THE  PRINI  LEVEL. 

now  wE  solve  the  problem  allowing  no  more 
than  200  STEPS,  ano  display  the  results  on 
the  PRINTER  (as  OPPOSED  TO  PLOTTING  THEM). 
NOTE  That  THE  LINE  EnO  FLAG  S  IS  NOT  RE¬ 
QUIRED  AFTER  THE  STOP  COMMAND. 


-The  same  job  in  more  compact  (and  less  readable)  form  appears  as  follows: 


HELP  1NIT  PRlN  10  FUNC  ROSE  POS I  ROSE  HESS  INIT  ANAL  UpOA  ANAL  BACK  S 
METH  PLSI  PARA  REAL  4,5000.  BACK  SEAR  CUBI  EVAL  PARA  DISP  BACK  BACK  S 
DISP  OPTI  BACK  ITER  1  INIT  PRIN  5  BACK  ITER  200  DlSP  RESU  PRIN  S 
STOP 
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;  Sec.  4.1  in  order  to  set  up  such  a  job.  However,  once  such  a  set  up  is 

made  for  a  particular  problem,  it  proves  very  convenient  for  executing  a 


number  of  studies. 

»  By  leaving  off  the  stop  and  insuring  that  the  last  line  does  not  end  in 


$,  a  packaged  command  sequence  similar  to  the  above  setup  is  also  very  helpful 
for  interactive  studies.  Simply  introduce  the  command  sequence  initially 
and  have  control  returned  to  an  interactive  terminal  afterwards  for  further 
processing. 

4.3  Examples 

In  this  section,  we  provide  sample  results  obtained  from  two  run 
streams.  The  first  is  taken  directly  from  an  interactive  session  at  a 
terminal.  The  commands  issued  by  the  user  are  underlined  for  emphasis. 


.JYOT  CA. 

HELLO*  1  AM  A  COMVEPImTIOMHL  PROGRAM  POP  LOCAL  FUNCTIONAL  MINIMIZATION. 
IF  YOU  MEED  HELP*  PLEA'.E  TYpE  HELP.  QTHERWI  SE  ... 

1  PL  EH  IE  T’i'FF  OmE  OF  THE  FOLLOWING  COMMAND:  ... 

INITIALIZE  ITERATE  DUPLAY 

pi  :p 

5  PLEA ' E  TYPE  ONE  OF  THE  FOLLOWING  COMMAND:  ... 

RESULTS:  DPTION'  COMMAND'  'TATE 

option: 


CUPPENT 

option: 

1 

POSITION 

RANDOM 

Zf 

k. 

FUNCTION 

SPECIAL 

.  3 

INITIAL 

HESSIAN 

ANAL  YT I C 

4 

UPDATE 

HESSIAN 

BF6S 

s 

PRI.NT 

e; 

•3 

METHOD 

PLSl 

i 

SEARCH 

QUADRATIC 

?  PLEASE  TYPE  ONE  OF  THE  FOLLOWING  COMMAND ‘  ... 
RESULTS  OPTIONS  COMMANDS  :TATE 

BACH  INIT 
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£  PLEASE 

TYPE 

ONE  DP  THE 

FOLLOWING  COMMANDS 

•  ■  • 

PRINT 

FUNCTION 

POSITION 

HF.SS  IhN 

METHOD 

eule 

SEARCH 

PARAMETER 

EVALUATE 

3  PLEASE 

TYPE 

ONE  OP  THE 

FOLLOWING  COMMANDS 

•  •  • 

RDSENBROC 

WOOD 

SPEC 

ELLIPSOID 

SPECIAL 

POWELL 

PLFASE  TYPE  TME  SIZE  <T>I  MEN?  I ON '  QF  YO'JP  FUNCT  I  ONAL  . 

4 

£  PLEASE  TYPE  ONE  OP  THE  FOLLOWING  COMMANDS  ... 

PRINT  FUNCTION  POSITION  HESSIAN 

METHOD  SEARCH  "'hPhMETEP'  EV'PLUnTE 

3  P-L EASE  TYPE  ONE  OF  THE  FOLLOWING  C CMMANP r  ... 
RDSENBROC  RANDOM  ZEPD  PERT' 

PPUELL  MOOD  RAYLEIGH 

RRYL  I  NGTDN 

INITIRL.  POSITION  -  PRYLEIGH 
1.000 00+ 00  4.71405-01  4.71405-01  4.71405-01 

£  PLER.'E  TYPE  ONE  OP  THE  FOLLOWING  COMMANDS  ... 

PRINT  FUNCTION  P'DIITIDN  HESSIAN 

METHOD  SEARCH  PARAMETER  E'.'RL  UATE 

EJCR_ 

FORM  OP  INITIRL  HE “SI RN  -  RNRLYTICRL 
1  PLER‘ E  TYpE  ONE  Or  THE  FOLLOWING  COMMAND ;  ... 

INITIRL IZE  1TEPRTE  DISPLAY 

ITEP  30 


♦  ♦♦  LDp£4  -  RN  UPDATE  COEFFICIENT  IS  ZEPO.  A*  F  =  0.  000 

PE-INITIRLIZE  THE  HESSIRN. 


POPM  OF  INITIRL  HESSIRN  -  RNRLYTICRL 
<G«  7  >  =  3. 7£91674-0£  RT  -4. 86871 1 £- 01  -  NDT 
fG-S-  =  -£. 0617467-04  RT  5. £478673-01  -  NDT 
1  PLEASE  TYPE  ONE  OP  THE  FOLLOWING  COMMANDS 
INITIALIZE  ITERRTE  DISPLAY 

PI~f  RESULTS 

3  PLEASE  TYPE  ONE  OP  THE  FOLLOWING  COMMANDS 
PPINT  PLOT 

PRINT 


CONVERGED 

CONVERGED 


o.  ooo 


PINAL  RESULTS 


FUNCTIONAL  VALUE 
GRADIENT  NORM 
STEP  NQPM 

ACCEPTANCE  CRITERION 
STEP  LENGTH  FACTOR 


1 . 0000000+00 
1 .6£?5?5£-07 
£. 3030643- 07 
3. £000000+01 
£. 0000000+00 
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NO. 

FUNCTIONAL  EVALS 

14 

NO. 

>?RAp]ENT  EVALS. 

34 

NO. 

HE  SSI AN  EVAl  S . 

d 

NO. 

HESSIAN  UPPhTES 

13 

NO. 

HES-VECTOP  MULTS 

3? 

NO. 

i tepations 

13 

NO. 

PE PEN.  PLANES 

0 

SOLUTION  (POSITION  QP  MINIMUM' 


1  #  f, 

356751-01  -3.6443101-03  3 

.6754551-03  -1 

.  3O03301-O3 

CONVERGENCE 

HISTORY 

rp 

cunstidnal 

pPAPI ENT  1 

CD :  *  THET 1 ^ »  1 ^ 

STEP 

ij 

i .  -•  ii  ij  n  u  1.14-  o  o 

1 . 715303+00 

1.303334-03 

1 . 3? 09 35+  0  0 

i 

1 . 664  al 1  +  00 

3.313740+00 

0.  000000 

0. 000000 

•• 

1 . 601733+00 

3.915311+00 

1.397474-03 

6. 741414-01 

o* 

1 . 314  34  3+00 

1 . 303314+00 

3.713347-07 

4.434947-01 

4 

1  .  041374  +  00 

6.  133373-01 

-7. 361933-06 

3.  346334-01 

C 

1  .  0  04503+00 

1.997154-01 

-4.704744-03 

3.  194359-01 

f 

1  .  0  00031 +00 

3.  347134-03 

-1 .403343-04 

4. 473335-03 

1 

1  .  i'i  0  0  i'i  03+00 

3.490150-03 

7.335336-05 

A.  93*' 069-03 

3 

1  .  0  00 0  0  0+ oo 

6.  9957  09-04 

3. 676903-03 

1 .413664-0  3 

—4 

1  .  0  0  0  0  o  0  +  o  o 

1 . 93964  0-04 

-3.517363-01 

1  .  163739-04 

1  0 

1  .  0  0 000  0+ oo 

3.433491-05 

4. 004605-03 

5. 3550"'?-05 

1 1 

1  .  0  00000+00 

4.753469-06 

3. 739343-04 

7.71 1701-06 

13 

1  .  0  0  0  0  0  0+  0  0 

1 .  167751-06 

-3. 4  0pb9.-i-03 

1 . 6777 36—  06 

1  3 

j .  oooo oo+ oo 

1 . 639595-07 

-5. 353430-03 

3. 303064-07 

P: 

_ea:e  type  dne 

□F  THE  FDLLDWING  COMMANDS 

•  •  • 

P& i nt  PLOT 

ip:-  option: 


CURRENT  OPTIONS 


1  POSITION 

RAYLEIGH 

3  ^UNCTION 

SPECIAL 

3  INITIAL 

HESSIAN 

ANALYTIC 

4  up p ATE 

5  PRINT 

HESSIAN 

5 

PL  SI 

PFAS 

6  METHOP 

7  SE APCH 

QUADRATIC 

PLEASE  TYPE 

□NE  OP  THE 

FOLLOWING  COMMANPS 

CULTS 

□PTIDNS 

commanps 

stop 


•  •  • 

STATE 


PUN  TcC'MINmTFIi 


The  second  example  is  the  result  of  submitting  the  run  stream  at  the 
beginning  of  Section  4.2  in  a  batch  processing  mode.  In  this  instance,  the 
user  commands  were  not  displayed  in  the  output.  Can  you  tell  where  they 
would  have  appeared  in  interactive  mode?  There  is  an  echo  option  for  batch 
mode  which  includes  a  display  of  the  commands. 


PXQT  CA. 

HELLO,  I  A*/!  A  CONVERSATIONAL  RR03SAV,  FOR  iCCAl  FUNCTIONAL  MINIMIZATION. 
IF  YOU  NEED  HELP.  PLEASE  TYPE  HELP.  OTHERWISE  ... 

1  PLEASE  TYPE  ONE  OF  THE  FOLLOWING  COMMANDS  ... 

INITIALIZE  iterate  display 

OUR  CONVERSATION  IS  ORGANIZED  AS  A  TREE  STRUCTURE.  AT  EACH  LEVEL  YOU 
WILL  HAVE  DIFFERENT  COMMAND  OPTIONS  TO  DICTATE  MY  SUBSEQUENT  ACTIONS. 
AT  EACH  LEVEL  I  SHALL  TElL  YOU  THE  CURRENT  LEVEL  NUMBER  AND  THE  VALID 

COMMANDS.  YOU  SElECT  ONE  AND  I  WILL  CARRY  IT  OUT.  YOU  MAY  MISSPELL  The 

COMMANDS  AS  LONG  AS  THE  FIRST  FOUR  CHARACTERS  ARE  CORRECT. 

WHEN  YOU  HAVE  COMPLETED  WORK  AT  A  BRANCH  OF  THE  TREE.  USE  THE  3ACKUP 

COMMAND  TO  GO  BACK  TO  THE  PREVIOUS  (HIGHER l  LEVEL  BRANCH.  FOR  BREVITY, 

I  DO  NOT  REPEATEDLY  PRINT  BACKUP  AS  A  COMMAND  OPTION.  THE  SAME  BREVITY 
CONVENTION  APPLIES  TO  THE  HELP  AND  STOP  COMMANDS. 


the  evaluation  phase  of  initialization  is  always  required  at  least  once 
BEFORE  The  SOLUTION  (ITERATION)  PHASE.  IT  SHOULD  BE  DONE  AFTER  MOST  OF 
THE  OTHER  INITIALIZATION  WORK  IS  DONE. 


GOOD  LUCK. 

INITIAL  POSITION  -  RDSENBROCK  STANDARD 

-1 .20000+00  1 . 00000*00 

GRADIENT 

-2.1560000*02  -B. 7999999*01 
F  «  2.4200000*01 

FORM  OF  INITIAL  HESSIAN  -  ANALYTICAL 

HESSIAN 

1 .330000+03 

4.000000*02  2.000000*02 

the  integer  parameters  are 


1 

1 

1 

1 

1  10 

1 

1 

1 

1 

2 

100 

5 

0 

0  1 

1 

1 

0 

0 

0 

0 

0 

0 

0  0 

0 

1 

3 

5 

9 

1  1 

13 

15 

17  1262 

2507 

3752 

4997 

THE  REAL  PARAMETERS  ARE 


3.200000+01 
0.000000 
1 .000000+06 
4.693321-01 
1 . 000000  +  06 
1.000000+06 


0.000000 
0.000000 
1 .000000+06 
1  .000000*06 
1  .000000*06 
1 .000000  +  06 


0.000000 
1 .280000+02 
2. 32B677+02 
1 . 000000+06 
1 .000000+06 


5.000000+03 
0.000000 
1 . 000000+06 
1 . 562050+00 
1 . 000000+06 


0 . 000000 
2.420000+01 
1 . 000000+06 
1 . 000000+06 
1 . 000000+06 
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CURRENT  OPTIONS 


1  POSITION 

R0SEN3R0C 

2  FUNCTION 

ROSENSROC 

3  INITIAL 

HESSIAN 

ANALYTIC 

4  UPDATE 

HESSIAN 

ANALYTIC 

5  PRINT 

10 

6  METHOD 

P  LS 1 

7  SEARCH 

CUBIC 

STEP 

2.4742007-02 

3.8061216-01 

F  =  4.731 7656+00 

GRADIENT 

-4.6415992+00 

-1 . 2383759-01 

POSITION 

-1 . 1752580+00 

1 . 3806122+00 

HESSIAN 

1  .107233+03 

4 . 701032+02 

2.000000+02 

final  RESULTS 

functional  value  0.0000000 

GRADIENT  NORM  0.0000000 

STEP  NORM  2 . 4607389-06 

ACCEPTANCE  CRITERION  3.2000000+01 

STEP  LENGTH  FACTOR  5.0000000+03 

NO.  FUNCTIONAL  EVALS  6 

NO.  GRADIENT  EVAlS.  8 

NO.  HESSIAN  EVALS.  E 

NO.  HESSIAN  UPDATES  0 

NO.  HES-VECTOR  IWULTS  14 

NO.  ITERATIONS  7 

NO.  DEGEN.  PLANES  0 


SOLUTION  (POSITION  OF  MINIMUM) 


1.0000000+00  t.ooooooo+oo 


CONVERGENCE  HI  STORY 


ITER  FUNCTIONAL  GRADIENT  COS  (  THF.  T  (  G ,  S  )  )  STEP 

0  2.420000+01  2 . 328677+02  4.693321-01  1.562050+00 

1  4.731786+00  4.6A3250+00  -9.146012-02  3.814155-01 

2  1.403741+03  1.363379+03  8.328175-01  4.943733+00 

3  5.727937-02  4.786486-01  -1.039726-04  3.747062+00 

4  3.281071-01  2.56'674+01  1.175821-01  4.357235-01 

5  1.554799-07  7.355650-03-4.906501-01  5.655543-02 

6  2.692291-12  5.475962-05-2.711999-02  7.997023-04 

7  0.000000  0.000000  0.000000  2.460739-06 

RUN  TERMINATED 
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